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Abstract 

In this paper we have modeled magnetic nanoparticle 
mathematically and have utilized a semi-classical approach 
based on Landau-Lifshitz Gilbert (LLG) equations. Our 
approach has based on Finite Element Method (FEM) as a 
numerical technique for finding approximate solutions of 
partial differential equations (PDE). The solution approach is 
based either on steady state problems, or rendering the PDE 
into an approximating system of ordinary differential 
equations, which are then numerically integrated using 
standard techniques. The proposed model has applied to 
three kind of magnetic nanoparticles including FluidMAG- 
DXS, FluidMAG-CMX and FluidMAG-PEI. Then Received 
Signal Strength (RSS) localization is utilized to find the true 
location of tumor where nanoparticles attached to it. The 
results of simulation prove that the error location can be less 
than 0.1 Cm. Researchers and clinicians can take advantage 
of these results in order for nanotechnology to make an 
impact in the diagnosis and treatment of malignancy. 
Furthermore the proposed model for magnetic nanoparticles 
can be used in improvement of drug delivery, cancer 
detection, diagnosis, imaging, and therapy while reducing 
toxicity. 
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Introduction 

Cancer is a major cause of death throughout the world, 
accounts for nearly one-quarter of total human 
mortality. Cancers treatment at an early stage of 
development allowing complete surgical removal is 
the most effective cancer therapy. Thus better method 
for detecting tumor at an early stage with low cost, 
noninvasive and precise techniques are the topic of 
extensive ongoing search and still is not fully achieved 
[1]. To this date, there are many imaging techniques 
that have been developed for cancer diagnosis, such as 
x-ray Computed Tomography (CT), Magnetic 
Resonance Imaging (MRI), and Positron Emission 
Tomography (PET). Although these imaging 



techniques are highly effective and widely utilized, 
they are expensive and their data interpretation is 
time-consuming. Furthermore in the cases that the 
patient is not allowed to expose in strong magnetic 
field, some of traditional approaches are not more 
practical. 

Cancer nanotechnology has the potential to 
significantly enhance current techniques to cancer 
detection, drug delivery, diagnosis, imaging, and 
therapy while reducing toxicity associated with 
traditional cancer therapy [2-3]. A wide variety of 
nanomaterials are currently under investigation and 
development for application relative to cancer 
nanotechnology including; polymers, dendrimers, 
lipids, organometallic, carbon based materials and 
magnetic nanoparticles [4-6]. In this case a variety of 
nanoscale magnetic Resonance Imaging (MR) contrast 
agents have been developed [7]. These particles have 
received significant attention, given the widespread 
use of MR in clinical practice currently [8]. 

In nanotechnology cancer treatment, one of the most 
hopeful techniques is the use of magnetic 
nanoparticles. Magnetic nanoparticles have diameters 
ranging from several nanometers to several hundred 
nanometers (nm). They also possess special properties 
useful in the design of new cancer detection systems 
[9]. Living cells are typically 10 pm in diameter and 
include of many smaller components like endoplasmic 
reticulum, mitochondrion and nucleus. Proteins are 
even smaller with a typical size of just 5 nm, which is 
comparable to the dimensions of the current smallest 
man-made nanoparticles [10]. The comparable size of 
nanoparticles and biological building blocks suggests 
the possibility of using nanoparticles as probes to spy 
on the working of intracellular machinery without 
introducing much interference [10]. Recently scientists 
have found that magnetic nanoparticles can be 
especially helpful in locating cancerous cell clusters 
during MRI scans. Like teeny guide missiles, the 
nanoparticles seek out tumor cells and attach 
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themselves to them. Once the nanoparticles bind 
themselves to these cancer cells, the particles operate 
like radio transmitters, greatly aiding the MRI's 
detection capability. Advancement in molecular 
biology is increasing understanding of rationale and 
relevant molecular targets in cancer. However, due to 
the heterogeneity of the molecular abnormalities in 
multiple tumor types, strategies designed to interfere 
with multiple molecular abnormalities will be 
necessary to impact survival. Nanoparticles have the 
potential to provide therapies not possible with other 
drug modalities. Researchers and clinicians must take 
advantage of these opportunities in order for 
nanotechnology to make an impact in the diagnosis 
and treatment of malignancy. 

In this paper we proposed a mathematical model for 
any kind of magnetic nanoparticles. Then by utilizing 
this model we can find the true location of tumor 
where nanoparticles are attaching themselves to it. 

This paper consists of four sections. In section II a 
mathematical model has been proposed. In section III 
the detector system is studied and Results of 
Simulations of our model are presented in section IV, 
which approves our theoretical work. 

Model of Magnetic Nanoparticles 

Problem Statement 

Spin and the related magnetic moment are primary 
properties of elementary particles. The spin quantum 
number relates to many phenomena in elementary 
particle physics as well as in nuclear and statistical 
physics. Spin is one type of angular momentum which 
does not have a classical analogy, but one may 
tentatively view a spin coming from a spinning 
motion of a nanoparticle. To this date, advancement in 
lithography and nanotechnology fabrication makes 
nanoparticles useful and feasible. In this case using of 
nanoparticles are needed a precise mathematical and 
physical model. To satisfy this demand we find a 
desirable model by utilizing Landau-Lifshitz Gilbert 
(LLG) equation. Then we perform this model to 
Received Signal Strength (RSS) localization system to 
find the location of nanoparticles. The received signals 
are coming from magnetic nanoparticles bound to 
"tumors". We aim to use this simulation study to 
reveal the theoretical possibility of magnetic 
nanoparticles for developing a low cost system 
capable of screening various types of cancers in a large 
population. 



In magnetic nanoparticles, quantum mechanical 
effects become significant and Maxwell's equations 
alone will not be adequate to describe the behavior of 
magnetic materials. Obviously, simulating these nano- 
magnetic materials from a quantum physics point of 
view would give the most perfect description. 

Due to have a mathematical model for magnetic 
nanoparticle we have utilized a semi-classical 
approach based on LLG equations [11]-[12] for nano 
scale materials. 

The Bloch equation describes well the magnetization 
dynamics [13] of an ensemble of non-interacting or 
weakly interacting spins, but it does not capture the 
proper physics of a strongly interacting spin system 
such as a piece of magnet. This is because the 
magnetization magnitude of a magnet shall not 
change with time, and does not preserve the 
magnitude of the magnetization. To take into account 
the dissipative effect of the environment, we utilize 
LLG equations to describe our system instead of Bloch 
equations. 

In our approach, we have used Finite Element Method 
(FEM) as a numerical technique for finding 
approximate solutions of partial differential equations 
(PDE). The solution approach is based either on steady 
state problems, or rendering the PDE into an 
approximating system of ordinary differential 
equations, which are then numerically integrated 
using standard techniques. 

Landau-Lifshitz Gilbert Equations 

Magnetization magnitude |M|for a specific material is 

constant and just dependant on environment 
temperature. The magnetization can be described as: 

M = \M.\-{mJ + m y ] + m z k) (1) 

Magnetic Gibbs free energy in nanoparticle (Etotal) 
consists of four elements: exchange energy, Zeemen 
energy, magnetostatic energy, and the anisotropy 
energy. These energies can be express as: 

E exeMn g e=^ m J + V m y? + V m z? dv ^ 

E magnetostatic= ^ H in- Mdv ^ 

E Zeemen =vjn ex -Mdv (4) 

Where, C is exchange stiffness, and for nanoparticles is 
in order of 10-13(J/m) [14]. H in is magnetic induction 



38 



Advances in Microelectronic Engineering (AIME) Volume 1 Issue 3, July 2013 



www .seipub .org/ aime 



for demagnetization field, H et is the magnetic 
induction vector of the externally applied field and p 
is the permeability of free space. Furthermore there is 
also the anisotropy energy, which is intrinsic to the 
material. The most common form of anisotropy is the 
stress anisotropy. Crystalline anisotropy occurs due to 
the spin-orbit interaction of the electrons. Depending 
on their lattice structure, the magnetic materials have a 
preferred direction of magnetization. This is called the 
easy axis. Anisotropy energy can be described by [15]: 

E anisotropy = L a ' V • sin 2 (/?) (5) 

Where (3 is the angle between the easy axis and the 
magnetization, La is the anisotropy constant, and V is 
the sample volume. 

The dynamic of magnetic properties is governed by 
LLG equation: 

(l + « 2 )— = -^VIxH r -«ymx(MxH r ) (6) 
dt 

Where y is the gyromagnetic ratio (in order of 10 5 ), a is 
the damping constant and H r is effective magnetic 
field: 

BF 



My — 



8M 



(7) 



Solving LLG equation for nanoparticles is required 2D 
mesh generator because nano magnetic materials 
currently under investigation are magnetic thin films. 

Thin films are really just a 2D structure where the 
thickness of the films is usually uniform. 

We use a mesh generation algorithm which is based 
on a mechanical correlation between a triangular mesh 
and a 2D truss structure. In this algorithm any set of 
points in the x, y, z-plane (nanoparticles are in z=0) can 
be mapped onto the "natural coordinates" (e, r\, E,) by 
the Delaunay algorithm [16]. In our proposed model 
the edges of Triangles (the connections between pairs 
of points) correspond to bars, and the points 
correspond to joints of the truss. 

Each bar has a force-displacement relationship f(L,L0) 
depending on its current length L and its unextended 
length Lo [16]. The magnitude of this force is just large 
enough to keep the node from moving outside. The 
positions of the joints (these positions are our principal 
unknowns) are found by solving for static force 
equilibrium in the desired structure. The hope is the 
lengths of all the bars at equilibrium will be nearly 
equal, giving a well-shaped triangular mesh. 

We can solve force equilibrium in our structure by 



considering all (x,y,z)-coordinates in mesh points into 
an N-by-2 array Q=[x y z]. So the force equilibrium 
equation could be written: 

F q (Q) = [FJQ) + F iy (Q) + F lz (Q}\+ [FJQ) + F ey (Q) + F ez (Q)\ (8) 

Where Fi is the internal forces from the bars, and Fe is 
the external force (reactions from the boundaries). The 
Delaunay algorithm find out non overlapping 
triangles that fill the convex hull of the input points, 
such that every edge is shared by at most two triangles, 
and the circumcircle of every triangle contains no 
other input points [16]. 

For solving force equilibrium equation, we consider 
the system of ODEs: 

dt ^ (9) 
Q(0 + ) = Q in 



The system is using forward Euler method for true 
approximation. The time is discretized by t n = nAt and 
solution Q n = Q(t n ) is update by: 

Q n+l =Q n +AtF{Q n ) (10) 

There are some force functions that have been studied 
in [17-18]. The function k(Lo,L) models ordinary linear 
springs[17]. Our implementation uses this linear 
response for the repulsive forces but it allows no 
attractive forces: 

\k{L -L) L>L Q 
I L < Ln 



f(L ,L) = - 



(11) 



If Lo is chosen constant then it generates uniform 
meshes. However there are many cases when it is 
advantageous to have different sizes in different 
regions. Where the geometry is more complex, it needs 
to be resolved by small elements (geometrical 
adaptively). The solution method may require small 
elements close to a singularity to give good global 
accuracy (adaptive solver). A uniform mesh with these 
small elements would require too many nodes. 

In our model we are using function H(x,y,z)as an 
element size to provide a desired edge length 
distribution. This function gives the relative 
distribution over the domain. So the number of node is 
not need to specify. In our implementation we use 
following function in MATLAB; function 
[Q,tri]=me_gen(G,L,h,B,fint,varargin). 

Q is N-by-3 array contains the(x, y, z) coordinates for 
each of the N nodes. Size of N is important parameters 
for our model. The triangle indices tri. The row 
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associated with each triangle has 3 integer entries to 
specify node numbers in that triangle. G is a distance 
function which is given for the geometry. G returns 
the signed distance from each node location Q to the 
closest boundary. 

The desired edge length function H(x, y, z) is given as 
a function L, which returns h(x,y,z) for all input points. 
The parameter h is the distance between points in the 
initial distribution Qo. For uniform meshes, the 
element size in the final mesh will usually be a little 
larger than this input [16]. 

The bounding box for the region is an array B=[xmin, 
ymin; xmax, ymax] (the particle is in Z-plate). The 
fixed node positions are given as an array pfix with 
two columns. Additional parameters to the functions 
G and L can be given in the last arguments varargin. 
The variable varargin is a cell array containing the 
optional arguments to the function. It must be 
declared as the last input argument and collects all the 
inputs from that point onwards. 

We will use directly with d = -j[x 2 + y 2 + z 2 ) - 1 which 

can be specified as an inline function. For a uniform 
mesh, h(x, y, z) returns a vector of l's. The circle has 
bounding box -1 < x < 1,-1 < y < 1, z = Z , with no fixed 

points. A mesh with element size approximately ho = 
0.2 is generated and the result is shown in Fig.l. 

Any variable cj), Including the x, y, and z coordinates, 
can then be interpolated within the proposed mesh 
grid by: 



i=\ 1 1 
And then Jacobian matrix for Eq. 12 is: 



(12) 



88 cW- 
Z 1 



88 cW- 



88 cW ; 



i=\ ds 1 i=\ ds 1 i=l ds 1 



88 cW- 
Z 1 



88 cW- 
Z 1 



88 cW- 

z ' 



i=\ drj 1 i=\ drj 1 i=\ drj 1 



88 cW- 
Z 1 



88 cW- 
Z 1 



88 cW- 

z ' 



(13) 



i=l d£ 1 i=\ d£ 1 i=\ d£ 1 

Applying Eq.12 to Eq.l gives: 

, . 88 88 88 

M = \M S | • ( S Nfm Yi i + Z Nfm vi j + Z N^k) (14) 



s\ v^.'T'm' ' i""yi- 
i=l i=l J 



i=\ 



Due to calculate the magnetic volume charge density 
we have: 



, V • \\M S • ( I N;mJ + S Nm :j + I Nm -k) 
1 . V ;=1 i=l - (=1 , 



\n V 



dv' (15) 



r - r 




FIG. 1 AN UNIFORM MESH WITH ELEMENT SIZE OF 0.2. FN 
THIS MESH THE SIZE OF P IS 88. 

And with some algebraic manipulation we have Eq. 16. 
Similarly, surface charge density is calculated by Eq. 
17. 



\M 



dN ; dN ; dN: 

(J— '-dv)m xj +(J — ldv)m ■ +(J— ldv)m d 
" " v dy ■ ' 



v dx 



v dz 



4,T 



'' ■ ! - ■ S dz 



(16) 



S dx 



S dy 



oo o 6o 9 bo o 



(17) 



For evaluating surface charge density, the unit normal 
vector should be calculated. Unit normal vector for 
any surface can be calculated by using Jacobian matrix. 

d(y, z) d(x, z) d(x, y) 

N(s,rj) = 1 + j + k (18) 

d(e, rj) d{s, rj) d(e, rj) 

The volume and surface charge terms are then 

summed for each element and then assembled into a 

global matrix just like any other finite element method: 

<P 'total =<PV +( P S (19) 

We have used interpolation functions as a weighting 
for Galerkin method. Then we multiply these 
weighting functions with LLG equation and integrate 
over the entire volume. 

We then obtain the formulation of the LLG equation 
by using Green's theorem, which states that: 
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dt 



i + a 
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M Lu ;=U=lv 



j=lk=l ™ jJ y J ^ y 



1+a MLu j'=U=1/=1v 



+m-m ,m ,+m-m ,m ,+m-m , m ,+2m-m ,m ,) 
x] zk zl xj yk yl z] xk zl xj xk xr 
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+ f =lk l^l (N i N i VN k ■ VN l dv}{m xj m z k m zl + m xj m xk m xl +m xj m yk m yl ]) 
88 88 88 

+ | 1 ^ 1/ J 1 Vxt (a v ( W* w »* ), ^- y v W^/* ) *)« ) 

88 88 88 
+^Z Z^m^ftJt^ 

88 88 88 

+ Z Z lm xj m k ((l(N i N j N k VN l dv)M-(lN i N j N l N,dv)h l ) 

j=\k=\l=\ J XK y * J K L v t J K I XI 

88 

-Z (((J (tf i W-*)-l)-(Jff.ff Mh ■)) 

j=l V J v 

Where; 

^ =$N i N j dv 

V 

In matrix term Eq. 20 is expressed as: 



(20) 



(21) 



(22) 



F(mx) is the derivative of rrix with respect to time Gi is 
the combined terms on the left hand side of Eq. 20. Wij 
is the coefficient for F(m x ) and it is a symmetrical 
matrix. We have been solved Eq. 22 by using time 
stepped integration and the fourth order Runge-Kutta 
formula [17]. In This method (m x , m y , m z ) should be 
normalized at the end of every time step iteration to 



guarantee: 



m x i +m j + m k 



1 



In Eq. 22 Wij is 88x88 matrix, and MATLAB cope 
with this kind of large matrix proficiently using the 
function sparse. Using the function sparse speed up 
the matrix inversion calculation and will make the 
whole system is a real time system. 

The global matrix for the FEM is assembled from these 
elemental matrices in the usual way. We use fourth- 
order Runge Kutta formula [18] to advance the 

solution from m to m 

n+l 

= m + 

6 3 3 6 

Where: 



n+l 



n ^1 ^2 *3 ^4 

m + — + — + — + - 



(23) 



fcj = (timestep)F(m ) 



M , 1 



&2 = (time step) F (m 

n ^2 

L = (timestep)F (m H ) 

3 2 

Yl 

k^ = (timestep)F(m +k^) 

In our simulation a time step on the order of 1 fs is 
usually provided a good accuracy. 

In Runge Kutta formula, the easiest way to stop the 
step time iteration is to terminate when the rate of 
change of the magnetization component is less than 
0.01° per fs. If the time step is in fs, then this rate of 
change is given by: 

n " n+1 180 

T )x <0.01 (24) 



d9 -[ m " ■ m 

— = cos ( 

dt 



n n+l 
m in 



n x timestep 



Generate niesh 
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Generate coefficient for scalar 
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Start time step 
integration 
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scalar potential 



Update magnetization 
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FIG. 2 THE FLEWCHART OF PROPOSED MODEL 
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The described algorithm has been illustrated in details 
in FIG. 2. 

2 

EJ 



EJ 



l 

EJ 

5 cm f 



3 

EJ 



4 

EJ 
EJ 

FIG. 3 THE DETECTOR SYSTEM INCLUDING 4 SQUIDs 



Detector System 

Our detection system is containing four gradient coils 
of SQUIDs (gradiometers) [19] marked by numbers 1 
through 4 (FIG. 3). 

The SQUIDs are used to detect and measure the 
signals which are produced by magnetic nanoparticles. 
They have a capability to detect vary small magnetic 
fields due to the high sensitivity [20]. 

Commercial Magneto-Encephalo-Graphic (MEG) 
systems, which are usually equipped with several 
hundreds of SQUID channels, are increasingly 
available in clinical settings. These systems are capable 
of measuring extremely weak magnetic fields on the 
order of 10-14T produced by active groups of neurons 
within the brain[20]. 

Applying the measured signals to the proposed 
algorithm four distances can be achieved. Here Di is 
indicated the distance between i-th coil and magnetic 
nanoparticles. The magnetic field signal from the 
tumor (magnetic nanoparticles) is calculated at the 
central point of the two coils. The spacing between 
each pair of the gradient coils is 5 cm wide and 0.75 
cm high. 

The detection system provides four distance [Di D2 D3 
D4]. Using RSS localization, the location of tumor can 
be identified precisely. 

Simulation and Result 

Our study utilized three types of magnetic 
nanoparticles manufactured by Chemicell GmbH in 
Germany as listed in table 1. 

The particles from Chemicell GmbH were supplied in 



the form of ferrofluid with dispersion of magnetic 
nanoparticles (iron oxides). These particles are coated 
with hydrophilic polymer matrices that prevent them 
from being aggregated by foreign ions. Using terminal 
functional groups of agents such as ion-exchange 
group or reactive group for covalent immobilization, 
these ferrofluids can bind to biological molecules as 
biomarkers for a number of applications such as MRI- 
based diagnosis and targeted drug delivery. 

TABLE 1 THREE KIND OF NANOPARTICLES USING IN THIS PAPER 



Product 
name 


Core 


Size 
(run) 


Polymer matrix 


FluidMAG- 
DXS 


FesCU 


100 


Dextran sulfate 


FluidMAG- 
CMX 


FesCU 


250 


Polyethylenimine 


FluidMAG- 
PEI 


FeiOi 


250 


Carboxymethyl dextran 



Magnetic field simulation: 

In this simulation we are applying the proposed 
algorithm which is illustrated in FIG. 2 over three 
magnetic nanoparticles. Then the simulation result of 
the algorithm is compared to our experimental data. 
The experimental data is achieved by using the 306- 
channel Elekta Neuromags MEG system 
manufactured by Elekta Neuromag Oy (Helsinki, 
Finland). 
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FIG. 4 SIMULATING FluidMAG-PEI WITH 1.212 DAMPING 
CONSTANT 
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FIG. 5 SIMULATING FluidMAG-CMX WITH 0.932 DAMPING 
CONSTANT 
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Tumor localization: 

In order to simulate the RSS localization based on the 
proposed algorithm we have utilized a simplified 
human model as shown in FIG. 7. 



* ex peri menial rc-sult 
proposed algorithm ri 



tt 



FIG. 6 SIMULATING FluidMAG-DXS WITH 0.932 DAMPING 
CONSTANT 



x (Cm) 




y(Cm ) 



FIG. 7 A SIMPLIFIED MODEL FOR HUMAN BODY 




Simulation, number 

FIG. 8 LOCALIZATION ERROR OVRE 1000 TEST USING 
FluidMAG-PEI 

Conclusion 

In this paper, we have studied a tumor localization 
system based on magnetic nanoparticles. The 
proposed system is consisting of four gradiometers. 
The SQUIDs which act as detectors sense the signals 
produced by the magnetic nanoparticles attached to 
tumors. We have proposed mathematical model for 
magnetic nanoparticles based on LLG equations. In 
this system nanoparticles act as a signal sources in the 
human body and four SQUIDs act as base station that 
find the true location of nanoparticles. 



Our simulation has provided estimates of the system 
sensitivity for tumors located at different sites of the 
body. The results of simulation prove that the location 
error in 1000 test is less than 0.1 Cm is around 0.08Cm. 

Our results have shown that, by combining 
nanotechnology with SQUIDs, tumors can be detected. 
Using this new method, it is possible to detect tumors 
at early stages containing only hundreds of cancerous 
cells. Once this system is developed, it may be used to 
scan a large number of people rapidly to detect 
various types of cancers at different body locations. 
The actual design and construction of this system as 
well as the tumor localization will be our next steps of 
investigation. 
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